library("easypackages")
packages("tidyverse","fpp3", "tsibble", "feasts","fable", "patchwork")
library(plotly)09 Resumen regresión dinámica
Modelos de regresión dinámica
1. Introducción
Los modelos de suavización exponencial y ARIMA no pueden incluir información externa a la serie. Por el otro lado, regresión lineal sí toma en cuenta otras variables, pero siempre asume que el error es ruido blanco. Por lo tanto, podemosunir estos dos modelos y resultará en un modelo de regresión dinámica, igual que al lineal pero el error es modelado con ARIMA.
Hay que tomar en cuenta que en regresión dinámica, es necesario que todas las variables sean estacionarias.
2. Regresión con errores ARIMA en R con fable
Al agregar pdq(d=1), R aplicará las primeras diferencias a todas las variables, si no se agrega los valores de pdq, R los buscará automáticamente.
ARIMA(y ~ x + pdq(1,1,0))<ARIMA model definition>Ejemplo
Analizar y pronosticar los cambios en el consumo personal a través de el ingreso disponible.
us_change <- read_csv("https://otexts.com/fpp3/extrafiles/us_change.csv") %>%
mutate(Time = yearquarter(Time)) %>%
as_tsibble(index = Time)Rows: 187 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): Consumption, Income, Production, Savings, Unemployment
date (1): Time
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.us_change %>%
gather("var", "value", Consumption, Income) %>%
ggplot(aes(x = Time, y = value)) +
geom_line() +
facet_grid(vars(var), scales = "free_y") +
xlab("Year") + ylab(NULL) +
ggtitle("Quarterly changes in US consumption and personal income")
us_change %>%
ggplot(aes(x = Income, y = Consumption))+
geom_point()
pdq automáticos:
fit <- us_change %>%
model(ARIMA(Consumption ~ Income))
report(fit)Series: Consumption
Model: LM w/ ARIMA(1,0,2) errors
Coefficients:
ar1 ma1 ma2 Income intercept
0.6922 -0.5758 0.1984 0.2028 0.5990
s.e. 0.1159 0.1301 0.0756 0.0461 0.0884
sigma^2 estimated as 0.3219: log likelihood=-156.95
AIC=325.91 AICc=326.37 BIC=345.29Viendo los errores de la regresión y de ARIMA con residuals(), los podemos seperar con type = “regression” y type = “innovations” para los errores ARIMA.
bind_rows(
`Regression Errors` = as_tibble(residuals(fit, type="regression")),
`ARIMA Errors` = as_tibble(residuals(fit, type="innovation")),
.id = "type"
) %>%
ggplot(aes(x = Time, y = .resid)) +
geom_line() +
facet_grid(vars(type), scales = "free_y") +
xlab("Year") + ylab(NULL)
Asegurarnos que los errores ARIMA sean ruido blanco:
fit %>% gg_tsresiduals()
Ljung-Box:
augment(fit) %>%
features(.resid, ljung_box, dof = 5, lag = 8)# A tibble: 1 × 3
.model lb_stat lb_pvalue
<chr> <dbl> <dbl>
1 ARIMA(Consumption ~ Income) 5.89 0.1173. Pronóstico
Para llevar a cabo pronósticos de modelos de regresión con errores ARIMA, se necesita realizar el pronóstico de:
- La parte de la regresión
- La parte de los erroes ARIMA
y combinar los resultados.
Necesitamos pronósticos de las variables independientes o predictoras para poder pronosticar nuestra variable de interés. Si las predictoras no son conocidas en el futuro, tenemos que modelarlas por separado, o asumir valores futuros para cada una.
Ejemplo:
Regresando al consumo personal a partir del ingreso, podemos suponer que los cambios porcentuales en el ingreso serán iguales a el cambio promedio porcentual de los últimos 40 años:
us_change_future <- new_data(us_change, 8) %>% mutate(Income = mean(us_change$Income))
forecast(fit, new_data = us_change_future) %>%
autoplot(slice_tail(us_change, n = 80)) +
xlab("Year") +
ylab("Percentage change") +
ggtitle("Pronóstico de regresión con errores ARIMA")
Comparando con ARIMA:
fit_prev <- us_change %>%
model(ARIMA(Consumption ~ PDQ(0,0,0)))
fit_prev %>% forecast(h=10) %>%
autoplot(slice(us_change, (n()-80):n())) +
ggtitle("Pronóstico con modelo ARIMA")
En este nuevo modelo, podemos capturar más informción, por lo tanto, los intervalos de predicción se reducen.
Hay que tomar en cuenta que los intervalos de predicción no toman en cuenta la incertidumbre de las predictoras.
Ejemplo:
Demanda de energía a parir de temperatura.
vic_elec_daily <- vic_elec %>%
filter(year(Time) == 2014) %>%
index_by(Date = date(Time)) %>%
summarise(
Demand = sum(Demand)/1e3,
Temperature = max(Temperature),
Holiday = any(Holiday)
) %>%
mutate(Day_Type = case_when(
Holiday ~ "Holiday",
wday(Date) %in% 2:6 ~ "Weekday",
TRUE ~ "Weekend"
))
vic_elec_daily %>%
ggplot(aes(x=Temperature, y=Demand, colour=Day_Type)) +
geom_point() +
ylab("Electricity demand (GW)") +
xlab("Maximum daily temperature")
Como vemos, es una función cuadratica y la energia demandada sí depende del día y la temperatura.
Gráficas de tiempo:
vic_elec_daily %>%
pivot_longer(cols = c(Demand, Temperature), names_to = "vars",
values_to = "value") %>%
ggplot(aes(x = Date, y = value)) +
geom_line() +
facet_wrap(~ vars, ncol = 1, strip.position = "right", scales = "free")
Al ser una relación cuadrática, ajustaremos el modelo cuadrático con errores ARIMA, y agregamos una variable para indicar si el día fue hábil o no.
fit <- vic_elec_daily %>%
model(ARIMA(Demand ~ Temperature + I(Temperature^2) + (Day_Type=="Weekday")))
report(fit)Series: Demand
Model: LM w/ ARIMA(2,1,2)(2,0,0)[7] errors
Coefficients:
ar1 ar2 ma1 ma2 sar1 sar2 Temperature
-0.1093 0.7226 -0.0182 -0.9381 0.1958 0.4175 -7.6135
s.e. 0.0779 0.0739 0.0494 0.0493 0.0525 0.0570 0.4482
I(Temperature^2) Day_Type == "Weekday"TRUE
0.1810 30.4040
s.e. 0.0085 1.3254
sigma^2 estimated as 44.91: log likelihood=-1206.11
AIC=2432.21 AICc=2432.84 BIC=2471.18fit %>% gg_tsresiduals()
augment(fit) %>%
features(.resid, ljung_box, dof = 8, lag = 14)# A tibble: 1 × 3
.model lb_stat lb_pv…¹
<chr> <dbl> <dbl>
1 "ARIMA(Demand ~ Temperature + I(Temperature^2) + (Day_Type ==… 28.4 7.91e-5
# … with abbreviated variable name ¹lb_pvalueAl momento de pronosticar, necesitamos datos respecto a la tempreatura futura y podemos hacer dos cosas:
- Conseguir los datos de algún pronóstico meteorológico y meterlos al modelo.
- Podemos pronosticar bajo un escenario. Asumiremos que la temperatura máxima se mantendrá constante a 26 grados.
vic_elec_future <- new_data(vic_elec_daily, 14) %>%
mutate(
Temperature = 26,
Holiday = c(TRUE, rep(FALSE, 13)),
Day_Type = case_when(
Holiday ~ "Holiday",
wday(Date) %in% 2:6 ~ "Weekday",
TRUE ~ "Weekend"
)
)
forecast(fit, vic_elec_future) %>%
autoplot(vic_elec_daily) + ylab("Electricity demand (GW)")
vic_elec_future <- new_data(vic_elec_daily, 14) %>%
mutate(
Temperature = c(seq(31,34),33,rep(34,3),rep(33,6)),
Holiday = c(TRUE, rep(FALSE, 13)),
Day_Type = case_when(
Holiday ~ "Holiday",
wday(Date) %in% 2:6 ~ "Weekday",
TRUE ~ "Weekend"
)
)
vic_elec_future# A tsibble: 14 x 4 [1D]
Date Temperature Holiday Day_Type
<date> <dbl> <lgl> <chr>
1 2015-01-01 31 TRUE Holiday
2 2015-01-02 32 FALSE Weekday
3 2015-01-03 33 FALSE Weekend
4 2015-01-04 34 FALSE Weekend
5 2015-01-05 33 FALSE Weekday
6 2015-01-06 34 FALSE Weekday
7 2015-01-07 34 FALSE Weekday
8 2015-01-08 34 FALSE Weekday
9 2015-01-09 33 FALSE Weekday
10 2015-01-10 33 FALSE Weekend
11 2015-01-11 33 FALSE Weekend
12 2015-01-12 33 FALSE Weekday
13 2015-01-13 33 FALSE Weekday
14 2015-01-14 33 FALSE Weekday forecast(fit, vic_elec_future) %>%
autoplot(vic_elec_daily) + ylab("Electricity demand (GW)")
4. Regresión armónica dinámica
Las series de Fourier es un método alternativo de los dummies muy bueno para periodos estacionales muy largos. Ya que ARIMA y ETS son métodos que estiman estacionalidad, pero pierde efectividad entre mayor sea el periodo estacional, y que, las dummies estiman estacionalidad de cualquier tamaño, pero si es un periodo muy grande, se necesitarán muchas variables, el uso de series de Fourier es buena opción.
Una serie de Fourier es una función básica trigonométrica, tiene un seno y un coseno. Al ser una función trigonométrica periódica, simula una estacionalidad y el uso de esta nos puede ayudar a pronosticar una serie estacional.
Como se comentó, una serie de Fourier se consituye de un seno y coseno, y se pueden utilizar tantos pares de senos y cosenos como se desea.
El número de pares es lo que cambiará la suavización o el ajuste del método a los datos de entrenamiento. Entre más pares, mejor ajuste (hay que tomar en cuenta el over-fitting).
El máximo absoluto de pares usados para un buen pronóstico es el periodo estacional \(m / 2\).
Ejemplos de estacionalidades largas:
- Datos diarios, donde pueden tener estacionalidad anual. (m = 365)
- Datos semanales (m = 52)
- Datos por cada media hora con estacionalidad diaria. (m = 48)
Ventajas:
- Permite cualquier tamaño de estacionalidad.
- Se pueden simular varios tipos de estacionalidad al mismo tiempo (con distintas frecuencias).
- La suavización del patrón estacional es controlado.
- La dinámica de corto plazo puede ser controlada a través de errores ARIMA.
Ejemplo
Gasto en comidas fuera de casa.
aus_cafe <- aus_retail %>%
filter(
Industry == "Cafes, restaurants and takeaway food services",
year(Month) %in% 2004:2018
) |> summarise(Turnover = sum(Turnover))
aus_cafe# A tsibble: 180 x 2 [1M]
Month Turnover
<mth> <dbl>
1 2004 Jan 1895.
2 2004 Feb 1766.
3 2004 Mar 1873.
4 2004 Apr 1874.
5 2004 May 1846
6 2004 Jun 1758.
7 2004 Jul 1878.
8 2004 Aug 1850.
9 2004 Sep 1922.
10 2004 Oct 1915.
# … with 170 more rowsautoplot(aus_cafe)Plot variable not specified, automatically selected `.vars = Turnover`
Probamos con varios niveles \(K\) (número de pares):
fit <- aus_cafe %>%
model(
`K = 1` = ARIMA(log(Turnover) ~ fourier(K = 1) + PDQ(0,0,0)),
`K = 2` = ARIMA(log(Turnover) ~ fourier(K = 2) + PDQ(0,0,0)),
`K = 3` = ARIMA(log(Turnover) ~ fourier(K = 3) + PDQ(0,0,0)),
`K = 4` = ARIMA(log(Turnover) ~ fourier(K = 4) + PDQ(0,0,0)),
`K = 5` = ARIMA(log(Turnover) ~ fourier(K = 5) + PDQ(0,0,0)),
`K = 6` = ARIMA(log(Turnover) ~ fourier(K = 6) + PDQ(0,0,0))
)
glance(fit)# A tibble: 6 × 8
.model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list>
1 K = 1 0.00175 317. -616. -615. -588. <cpl [2]> <cpl [3]>
2 K = 2 0.00107 362. -700. -698. -661. <cpl [5]> <cpl [1]>
3 K = 3 0.000761 394. -763. -761. -725. <cpl [3]> <cpl [1]>
4 K = 4 0.000539 427. -822. -818. -771. <cpl [1]> <cpl [5]>
5 K = 5 0.000317 474. -919. -917. -875. <cpl [2]> <cpl [0]>
6 K = 6 0.000316 474. -920. -918. -875. <cpl [0]> <cpl [1]>rep_model <- function(modelo) {
print(strrep("- - ",18))
print(paste("Modelo", modelo, sep = " "))
fit %>%
select(all_of(modelo)) %>%
report()
}
names(fit) %>%
map(rep_model)[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 1"
Series: Turnover
Model: LM w/ ARIMA(2,1,3) errors
Transformation: log(Turnover)
Coefficients:
ar1 ar2 ma1 ma2 ma3 fourier(K = 1)C1_12
-0.2266 0.0471 -0.6371 -0.5810 0.6407 0.0067
s.e. 0.1119 0.1103 0.0817 0.0682 0.0651 0.0018
fourier(K = 1)S1_12 intercept
-0.0358 0.0041
s.e. 0.0018 0.0011
sigma^2 estimated as 0.001747: log likelihood=317.24
AIC=-616.47 AICc=-615.41 BIC=-587.78
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 2"
Series: Turnover
Model: LM w/ ARIMA(5,1,1) errors
Transformation: log(Turnover)
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 fourier(K = 2)C1_12
-0.7492 -0.6512 -0.0657 0.2450 0.5088 -0.2694 0.0067
s.e. 0.1299 0.1497 0.1650 0.1223 0.0750 0.1442 0.0020
fourier(K = 2)S1_12 fourier(K = 2)C2_12 fourier(K = 2)S2_12 intercept
-0.0363 -0.0044 -0.0217 0.0039
s.e. 0.0020 0.0015 0.0015 0.0010
sigma^2 estimated as 0.001073: log likelihood=361.85
AIC=-699.71 AICc=-697.83 BIC=-661.46
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 3"
Series: Turnover
Model: LM w/ ARIMA(3,1,1) errors
Transformation: log(Turnover)
Coefficients:
ar1 ar2 ar3 ma1 fourier(K = 3)C1_12
-0.4021 0.2096 0.4883 -0.6694 0.0067
s.e. 0.2453 0.2633 0.1475 0.2531 0.0022
fourier(K = 3)S1_12 fourier(K = 3)C2_12 fourier(K = 3)S2_12
-0.0363 -0.0045 -0.0215
s.e. 0.0022 0.0014 0.0014
fourier(K = 3)C3_12 fourier(K = 3)S3_12 intercept
-0.0071 -0.0355 0.0039
s.e. 0.0016 0.0016 0.0009
sigma^2 estimated as 0.0007609: log likelihood=393.61
AIC=-763.21 AICc=-761.33 BIC=-724.96
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 4"
Series: Turnover
Model: LM w/ ARIMA(1,1,5) errors
Transformation: log(Turnover)
Coefficients:
ar1 ma1 ma2 ma3 ma4 ma5 fourier(K = 4)C1_12
-0.7133 -0.0944 -0.0834 0.1182 -0.2059 0.3678 0.0067
s.e. 0.0806 0.1011 0.0720 0.0927 0.0642 0.0671 0.0018
fourier(K = 4)S1_12 fourier(K = 4)C2_12 fourier(K = 4)S2_12
-0.0364 -0.0043 -0.0215
s.e. 0.0018 0.0019 0.0019
fourier(K = 4)C3_12 fourier(K = 4)S3_12 fourier(K = 4)C4_12
-0.0069 -0.0355 0.0034
s.e. 0.0012 0.0012 0.0019
fourier(K = 4)S4_12 intercept
-0.0203 0.0039
s.e. 0.0019 0.0011
sigma^2 estimated as 0.0005386: log likelihood=426.78
AIC=-821.57 AICc=-818.21 BIC=-770.57
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 5"
Series: Turnover
Model: LM w/ ARIMA(2,1,0) errors
Transformation: log(Turnover)
Coefficients:
ar1 ar2 fourier(K = 5)C1_12 fourier(K = 5)S1_12
-0.4448 -0.1193 0.0068 -0.0364
s.e. 0.0741 0.0747 0.0024 0.0024
fourier(K = 5)C2_12 fourier(K = 5)S2_12 fourier(K = 5)C3_12
-0.0044 -0.0216 -0.0070
s.e. 0.0014 0.0014 0.0013
fourier(K = 5)S3_12 fourier(K = 5)C4_12 fourier(K = 5)S4_12
-0.0355 0.0035 -0.0202
s.e. 0.0013 0.0014 0.0014
fourier(K = 5)C5_12 fourier(K = 5)S5_12 intercept
-0.0069 -0.0249 0.0039
s.e. 0.0014 0.0014 0.0008
sigma^2 estimated as 0.0003173: log likelihood=473.73
AIC=-919.47 AICc=-916.91 BIC=-874.85
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 6"
Series: Turnover
Model: LM w/ ARIMA(0,1,1) errors
Transformation: log(Turnover)
Coefficients:
ma1 fourier(K = 6)C1_12 fourier(K = 6)S1_12 fourier(K = 6)C2_12
-0.3954 0.0068 -0.0363 -0.0044
s.e. 0.0619 0.0024 0.0024 0.0016
fourier(K = 6)S2_12 fourier(K = 6)C3_12 fourier(K = 6)S3_12
-0.0215 -0.0070 -0.0355
s.e. 0.0016 0.0014 0.0014
fourier(K = 6)C4_12 fourier(K = 6)S4_12 fourier(K = 6)C5_12
0.0035 -0.0202 -0.0069
s.e. 0.0013 0.0013 0.0013
fourier(K = 6)S5_12 fourier(K = 6)C6_12 intercept
-0.0249 0.0014 0.0039
s.e. 0.0013 0.0009 0.0008
sigma^2 estimated as 0.0003163: log likelihood=474.03
AIC=-920.06 AICc=-917.5 BIC=-875.44[[1]]
# A mable: 1 x 1
`K = 1`
<model>
1 <LM w/ ARIMA(2,1,3) errors>
[[2]]
# A mable: 1 x 1
`K = 2`
<model>
1 <LM w/ ARIMA(5,1,1) errors>
[[3]]
# A mable: 1 x 1
`K = 3`
<model>
1 <LM w/ ARIMA(3,1,1) errors>
[[4]]
# A mable: 1 x 1
`K = 4`
<model>
1 <LM w/ ARIMA(1,1,5) errors>
[[5]]
# A mable: 1 x 1
`K = 5`
<model>
1 <LM w/ ARIMA(2,1,0) errors>
[[6]]
# A mable: 1 x 1
`K = 6`
<model>
1 <LM w/ ARIMA(0,1,1) errors>p <- fit %>%
forecast(h = "2 years") %>%
autoplot(aus_cafe) +
facet_wrap(vars(.model), ncol = 2) +
guides(colour = FALSE) +
geom_label(
aes(x = yearmonth("2007 Jan"), y = 4250, label = paste0("AICc = ", format(AICc))),
data = glance(fit)
)Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
of ggplot2 3.3.4.p
p <- fit %>%
select(`K = 3`) %>%
forecast(h = "2 years") %>%
autoplot(aus_cafe) +
facet_wrap(vars(.model), ncol = 2) +
guides(colour = FALSE)
p
Ejemplo 2
vic_elec# A tsibble: 52,608 x 5 [30m] <Australia/Melbourne>
Time Demand Temperature Date Holiday
<dttm> <dbl> <dbl> <date> <lgl>
1 2012-01-01 00:00:00 4383. 21.4 2012-01-01 TRUE
2 2012-01-01 00:30:00 4263. 21.0 2012-01-01 TRUE
3 2012-01-01 01:00:00 4049. 20.7 2012-01-01 TRUE
4 2012-01-01 01:30:00 3878. 20.6 2012-01-01 TRUE
5 2012-01-01 02:00:00 4036. 20.4 2012-01-01 TRUE
6 2012-01-01 02:30:00 3866. 20.2 2012-01-01 TRUE
7 2012-01-01 03:00:00 3694. 20.1 2012-01-01 TRUE
8 2012-01-01 03:30:00 3562. 19.6 2012-01-01 TRUE
9 2012-01-01 04:00:00 3433. 19.1 2012-01-01 TRUE
10 2012-01-01 04:30:00 3359. 19.0 2012-01-01 TRUE
# … with 52,598 more rowsUsando datos cada hora:
(vic_elec_hour <- index_by(vic_elec, time = ~ floor_date(., unit="hour")) |>
summarise(demand = sum(Demand),
mean_temp = mean(Temperature),
max_temp = max(Temperature),
holiday = any(Holiday)) |>
mutate(
day_type = case_when(
holiday ~ "holiday",
wday(time) %in% 2:6 ~ "weekday",
TRUE ~ "weekend"
)))# A tsibble: 26,304 x 6 [1h] <Australia/Melbourne>
time demand mean_temp max_temp holiday day_type
<dttm> <dbl> <dbl> <dbl> <lgl> <chr>
1 2012-01-01 00:00:00 8646. 21.2 21.4 TRUE holiday
2 2012-01-01 01:00:00 7927. 20.6 20.7 TRUE holiday
3 2012-01-01 02:00:00 7902. 20.3 20.4 TRUE holiday
4 2012-01-01 03:00:00 7256. 19.8 20.1 TRUE holiday
5 2012-01-01 04:00:00 6793. 19.0 19.1 TRUE holiday
6 2012-01-01 05:00:00 6636. 18.7 18.8 TRUE holiday
7 2012-01-01 06:00:00 6548. 18.7 18.8 TRUE holiday
8 2012-01-01 07:00:00 6865. 19.6 20.1 TRUE holiday
9 2012-01-01 08:00:00 7300. 21.8 22.6 TRUE holiday
10 2012-01-01 09:00:00 8002. 24.6 25.2 TRUE holiday
# … with 26,294 more rowsp <- autoplot(vic_elec_hour, demand, color = "blue")
ggplotly(p, dynamicTicks = TRUE) |> rangeslider()vic_elec_hour |>
ggplot(aes(x = mean_temp, y = demand,
color = holiday)) +
geom_point(alpha = 0.5)
vic_elec_hour |>
ggplot(aes(x = mean_temp, y = demand,
color = day_type)) +
geom_point(alpha = 0.5)
(train <- filter(vic_elec_hour, year(time) %in% 2012:2013))# A tsibble: 17,544 x 6 [1h] <Australia/Melbourne>
time demand mean_temp max_temp holiday day_type
<dttm> <dbl> <dbl> <dbl> <lgl> <chr>
1 2012-01-01 00:00:00 8646. 21.2 21.4 TRUE holiday
2 2012-01-01 01:00:00 7927. 20.6 20.7 TRUE holiday
3 2012-01-01 02:00:00 7902. 20.3 20.4 TRUE holiday
4 2012-01-01 03:00:00 7256. 19.8 20.1 TRUE holiday
5 2012-01-01 04:00:00 6793. 19.0 19.1 TRUE holiday
6 2012-01-01 05:00:00 6636. 18.7 18.8 TRUE holiday
7 2012-01-01 06:00:00 6548. 18.7 18.8 TRUE holiday
8 2012-01-01 07:00:00 6865. 19.6 20.1 TRUE holiday
9 2012-01-01 08:00:00 7300. 21.8 22.6 TRUE holiday
10 2012-01-01 09:00:00 8002. 24.6 25.2 TRUE holiday
# … with 17,534 more rows# maximos:
# year = 4380
# week = 84
# day = 12
fit <- train |>
model(
Fourier = TSLM(log(demand) ~ fourier(period = "year", K = 438) +
fourier(period = "week", K = 42) +
fourier(period = "day", K = 6)))
glance(fit)# A tibble: 1 × 15
.model r_squa…¹ adj_r…² sigma2 stati…³ p_value df log_lik AIC AICc
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
1 Fourier 0.876 0.869 0.00444 122. 0 961 23128. -94120. -94008.
# … with 5 more variables: BIC <dbl>, CV <dbl>, deviance <dbl>,
# df.residual <int>, rank <int>, and abbreviated variable names ¹r_squared,
# ²adj_r_squared, ³statisticfc <- forecast(fit, h = "1 year")Warning: There was 1 warning in `mutate()`.
ℹ In argument: `Fourier = (function (object, ...) ...`.
Caused by warning:
! prediction from a rank-deficient fit may be misleadingfc# A fable: 8,766 x 4 [1h] <Australia/Melbourne>
# Key: .model [1]
.model time demand .mean
<chr> <dttm> <dist> <dbl>
1 Fourier 2014-01-01 00:00:00 t(N(9, 0.0047)) 8399.
2 Fourier 2014-01-01 01:00:00 t(N(9, 0.0047)) 8329.
3 Fourier 2014-01-01 02:00:00 t(N(9, 0.0047)) 8063.
4 Fourier 2014-01-01 03:00:00 t(N(8.9, 0.0047)) 7697.
5 Fourier 2014-01-01 04:00:00 t(N(8.9, 0.0047)) 7484.
6 Fourier 2014-01-01 05:00:00 t(N(9, 0.0047)) 7727.
7 Fourier 2014-01-01 06:00:00 t(N(9.1, 0.0047)) 8572.
8 Fourier 2014-01-01 07:00:00 t(N(9.2, 0.0047)) 9740.
9 Fourier 2014-01-01 08:00:00 t(N(9.3, 0.0047)) 10537.
10 Fourier 2014-01-01 09:00:00 t(N(9.3, 0.0047)) 10569.
# … with 8,756 more rowsggplot() +
geom_line(vic_elec_hour, mapping = aes(x = time, y = demand)) +
geom_line(fc, mapping = aes(x = time, y = .mean), color = "blue", alpha = 0.8)